home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1995 April / macformat-023.iso / Shareware City / Developers / RLaB / toolbox / ode4.r < prev    next >
Encoding:
Text File  |  1994-09-20  |  2.5 KB  |  140 lines  |  [TEXT/RLAB]
  1. //
  2. //  This is an adaptive 4th-order Runge-Kutta integrator.
  3. //  See Numerical Recipes.
  4. //
  5.  
  6. static (rk4, collapse);
  7.  
  8. ode4 = function( ydot, t0, tf, y0, output, tol, h, yscale ) 
  9. {
  10.   global (eps)
  11.  
  12.   TINY = 1e-4;
  13.  
  14.   //  Supply defaults.
  15.   if (!exist (tol)) { tol = 1e-3; }
  16.   if (!exist (h)) { h = abs( tf - t0 ) / 1000.0; }
  17.  
  18.   //  Initialize.
  19.   t = t0;
  20.   h = h * ( tf - t0 ) ./ abs( tf - t0 );
  21.   y = y0[:];
  22.   I = 2;
  23.  
  24.   // Save the initial  conditions
  25.  
  26.   if (!exist (output)) 
  27.   {
  28.     yout.[1] = [t, y.'];
  29.   else
  30.     yout.[1] = output(t,y).';
  31.   }
  32.  
  33.   while ( (t-tf)*(tf-t0) < 0.0 ) 
  34.   {
  35.  
  36.     dydt = ydot( t, y );
  37.     if (!exist (yscale)) { scale = abs(y) + abs(h*dydt) + TINY; }
  38.  
  39.     //    If the step can overshoot end, cut the stepsize.
  40.  
  41.     if ( (t+h-tf)*(t+h-t0) > 0.0 ) { h = tf - t; }
  42.  
  43.     //    Take a step
  44.     ys = y;
  45.  
  46.     while ( 1 ) 
  47.     {
  48.       //  Take two half steps.
  49.  
  50.       hh = 0.5 * h;
  51.       yt = rk4( ydot, t, hh, ys, dydt );
  52.       y =  rk4( ydot, t+hh, hh, yt, ydot( t+hh, yt ) );
  53.  
  54.       //  Check for step too small.
  55.  
  56.       if ( t+h == t ) 
  57.       {
  58.           fprintf("stderr", "rlab: run time error: Step size too small in ode4\n");
  59.           yt = 0.0;
  60.         errmax = 1.0;
  61.         tf = t;
  62.       else
  63.     //  Take the full step.
  64.           yt = y - rk4( ydot, t, h, ys, dydt );
  65.     
  66.     //  Evaluate accuracy.
  67.         errmax = norm( abs( yt ) ./ scale, "I" ) / tol;
  68.       }
  69.  
  70.       //  Success?
  71.       if ( errmax <= 1.0 ) { break; }
  72.  
  73.       //  Try again.  Reduce step size, but by no more than
  74.       //  a factor of about 50.
  75.  
  76.       if ( errmax > 4e6 ) 
  77.       {
  78.           h = 0.02*h;
  79.       else
  80.         h = 0.9*h*errmax.^-0.25;
  81.       }
  82.     }
  83.  
  84.     t = t + h;
  85.  
  86.     //    Take care of 5th order.
  87.     y = y + yt./15.0;
  88.  
  89.     //    Save output at the current point.
  90.  
  91.     if (!exist (output))
  92.     {
  93.       yout.[I] = [t, y.'];
  94.     else
  95.       yout.[I] = output( t, y ).';
  96.     }
  97.     I++;
  98.  
  99.     //    Determine next step size.
  100.     if ( errmax > 6e-4 ) 
  101.     {
  102.       h = 0.9*h*errmax.^-0.2;
  103.     else
  104.       h = 4*h;
  105.     }
  106.   }
  107.   return collapse (yout);
  108. };
  109.  
  110. rk4 = function(ydot, t, h, y, dydt) 
  111. {
  112.   local( k1, k2, k3, k4 );
  113.  
  114.   k1 = h * dydt;
  115.   k2 = h * ydot( t+h/2, y+k1/2 );
  116.   k3 = h * ydot( t+h/2, y+k2/2 );
  117.   k4 = h * ydot( t+h, y+k3 );
  118.  
  119.   return y + ( k1 + 2.0*k2 + 2.0*k3 + k4 ) / 6.0;
  120. };
  121.  
  122. //
  123. // Collapse a LIST of matrices into a single matrix.
  124. //
  125.  
  126. collapse = function ( LIST )
  127. {
  128.   local (i, m, row, col);
  129.  
  130.   row = size (LIST.[members (LIST)[1]])[1];
  131.   col = size (LIST.[members (LIST)[1]])[2];
  132.   m = zeros (length (LIST)*row, col);
  133.  
  134.   for (i in 1:length (LIST))
  135.   {
  136.     m[i;] = LIST.[i];
  137.   }
  138.   return m;
  139. };
  140.